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The  diagonal  forms  are  constructed  for  the  translation  operators  for  the  Helmholtz  equation 
in  three  dimensions.  While  the  operators  themselves  have  a  fairly  complicated  structure 
(described  somewhat  incompletely  by  the  classical  addition  theorems  for  the  Bessel  func¬ 
tions),  their  diagonal  forms  turn  out  to  be  quite  simple.  These  diagonal  forms  are  realized 
as  generalized  integrals,  possess  straightforward  physical  interpretations,  and  admit  stable 
numerical  implementation.  This  paper  uses  the  obtained  analytical  apparatus  to  construct 
an  algorithm  for  the  rapid  application  to  arbitrary  vectors  of  matrices  resulting  from  the 
discretization  of  integral  equations  of  the  potential  theory  for  the  Helmholtz  equation  in 
three  dimensions.  It  is  an  extension  to  the  three-dimensional  case  of  the  results  of  [13], 
where  a  similar  apparatus  is  developed  in  the  two-dimensional  case. 
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Diagonal  Forms  of  Translation  Operators  for  the  Helmholtz  Equation 

in  Three  Dimensions 


1.  Introduction 

One  of  standard  approaches  to  numerical  treatment  of  boundary  value  problems  for  ellip¬ 
tic  partial  differential  equations  (PDEs)  calls  for  converting  them  into  second  kind  integral 
equations  (SKIEs)  with  subsequent  discretization  of  the  latter  via  appropriate  quadrature  for¬ 
mulae.  Discretization  of  the  resulting  SKIEs  usually  leads  to  dense  large-scale  systems  of 
linear  algebraic  equations,  which  are  in  turn  solved  by  means  of  some  iterative  technique,  such 
as  Generalized  Conjugate  Residual  algorithm.  Most  iterative  schemes  for  the  solution  of  linear 
systems  of  this  type  require  application  of  the  matrix  of  the  system  to  a  sequence  of  recursively 
generated  vectors.  Applying  a  dense  matrix  to  a  vector  is  am  order  n2  procedure,  where  n  is  the 
dimension  of  the  matrix,  which  in  this  case  is  equal  to  the  number  of  nodes  in  the  discretization 
of  the  domain  of  the  integral  equation.  As  a  result,  the  whole  process  is  at  least  of  the  order 
n2,  and  for  many  large  scale  problems,  this  estimate  is  prohibitively  large. 

During  the  last  several  years,  a  group  of  algorithms  has  been  introduced  for  the  rapid  ap¬ 
plication  to  arbitrary  vectors  of  matrices  resulting  from  the  discretization  of  integral  equations 
from  several  areas  of  applied  mathematics.  The  schemes  include  the  Fast  Multipole  Method  for 
the  Laplace  equation  in  two  and  three  dimensions  (see,  for  example,  [7]),  the  fast  Gauss  trans¬ 
form  (see  [9]),  the  Fast  Laplace  Transform  (see  [13, 16]),  and  several  other  schemes.  In  all  cases, 
the  resulting  algorithms  have  asymptotic  CPU  time  estimates  of  either  0(n)  or  0(n  •  log(n)), 
and  are  a  dramatic  improvement  over  the  classical  ones  for  large-scale  problems. 

All  such  schemes  are  based  on  one  of  two  approaches. 

1.  The  first  approach  utilizes  the  fact  that  the  kernel  of  the  integral  operator  to  be  applied 
is  smooth  (away  from  the  diagonal  or  some  other  small  part  of  the  matrix),  and  decomposes  it 
into  some  appropriately  chosen  set  of  functions  (Chebychev  polynomials  in  [13]  and  [2],  wavelets 
in  [4],  wavelet-like  objects  in  [3],  etc.).  This  approach  is  extremely  general  and  easy  to  use, 
since  a  single  scheme  is  applicable  to  a  wide  class  of  operators. 
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2.  The  second  approach  is  restricted  to  the  cases  when  the  integral  operator  has  some  special 
analytical  structure,  and  uses  the  corresponding  special  functions  (multipole  expansions  for  the 
Laplace  equation  in  [7],  [8],  Hermite  polynomials  in  [9],  Laguerre  polynomials  in  [16],  etc. 

In  this  approach,  a  special-purpose  algorithm  has  to  be  constructed  for  each  narrow  class  of 
kernels,  and  in  each  case  the  appropriate  special  functions  and  translation  operators  for  them 
(historically  known  as  Addition  Theorems)  have  to  be  available.  However,  once  constructed, 
such  algorithms  tend  to  be  extremely  efficient.  In  addition,  there  are  several  important  sit¬ 
uations  where  the  first  approach  fails,  but  the  second  can  be  used  (a  typical  example  is  the 
n-body  graviational  problem  with  a  highly  non-uniform  distribution  of  particles,  as  in  [5]). 

Both  of  the  above  approaches  fail  when  the  kernel  is  highly  oscillatory,  and  simple  counter¬ 
examples  show  that  it  is  impossible  to  construct  a  scheme  that  would  work  in  the  general 
oscillatory  case  (the  Nyquist  theorem  being  the  basic  obstacle).  However,  several  oscillatory 
problems  are  of  sufficient  importance  that  it  is  worth-while  to  construct  special  purpose  meth¬ 
ods  for  them.  A  typical  example  are  kernels  satisfying  the  Helmholtz  equation  in  two  and  three 
dimensions,  since  this  is  the  equation  controlling  the  propagation  of  acoustic  and  electromag¬ 
netic  waves,  and  many  quantum-mechanical  phenomena.  Unlike  the  non-oscillatory  case,  the 
oscillatory  one  requires  a  fairly  subtle  mathematical  apparatus,  and  for  the  Helmholtz  equation 
in  two  dimensions,  such  an  apparatus  is  constructed  in  [14]. 

The  present  paper  presents  an  extension  of  the  results  of  [14]  to  the  three-dimensional  case, 
and  a  description  of  an  algorithm  for  the  rapid  application  to  arbitrary  vectors  of  matrices 
resulting  from  the  discretization  of  integral  equations  of  the  potential  theory  for  the  Helmholtz 
equation  in  three  dimensions.  The  principal  purpose  of  this  paper  are  Theorems  3.1  -  3.3,  de¬ 
scribing  the  diagonal  forms  of  the  well-known  translation  operators  for  the  Helmholtz  equation 
in  three  dimensions. 

2.  Analytical  and  Numerical  Preliminaries 

2.1.  Notation.  We  will  be  denoting  by  (p,9,<f>)  the  spherical  coordinates  in  R 3,  with  the 
Euclidean  coordinates  denoted  by  (x,y,  z).  Given  a  point  s  on  the  two-dimensional  sphere  52, 
we  will  denote  its  spherical  coordinates  by  ( 0(s),4>($ )),  and  note  that  the  north  pole  has 
the  coordinates  (jr,*),  while  the  coordinates  of  the  south  pole  s$  are  (0,*). 
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We  will  denote  by  E  the  natural  embedding  S 2  — ►  R3,  defined  by  the  formula 


E(s)  —  (cos(0(s))  •  cos(^(s)),co.s(0(.s))  •  stn(<£(s)),.stn(0(s))).  (1) 

For  a  non-zero  vector  u  €  R3,  we  will  denote  by  P(u)  the  point  on  S2  defined  by  the  formula 


u 


P(”)  - 


(2) 


Sometimes,  we  will  use  a  more  invariant  notation,  saying  that  the  pair  (r  €  Rl,s  G  52)  is  the 
spherical  coordinates  of  the  point  u  €  R3,  with  r,  s  defined  by  the  formulae 


r  =  u 


s  =  P(«). 


(3) 


(4) 


For  a  pair  of  points  sq,s  €  S2,  we  will  denote  by  (s)  the  angle  between  the  vectors  E(sq), 
E(s). 

Finally,  for  any  so>  *  €  S2,  we  will  denote  by  c(so,  s)  the  cosine  of  the  angle  between  the 
vectors  E(so),  E(s),  so  that 

c(s0,  s)  =  (E(s0),  E(s))  =  cos(Aao(s)).  (5) 


2.2.  Charges  and  dipoles.  For  a  Helmholtz  equation 

V2/  +  k2f  =  0  (6) 

we  will  define  the  potential  /*0  :  R3  \  {x0}  — ►  C1  of  a  unit  charge  located  at  the  point  xo  €  R3 
by  the  formula 

&(*)  =  M»ll*  -*.11),  (7) 

where  ho  denotes  the  spherical  Hankel  function  of  order  zero  (see  (16)  below).  For  any  h  €  R3 
such  that  ||h||  =  1,  we  will  define  the  potential  /* h  of  a  unity  dipole  located  at  xq  and  oriented 
in  the  direction  h  by  the  formula 

/£,»(*)  =  -'■.(*11*  -  *»li)  ■  (8) 
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As  is  well  known,  both  potentials  /*0,  /*0>A  (as  well  as  most  other  physically  meaningful 
potentials  for  the  equation  (6)),  satisfy  the  radiation  condition  at  oo,  i.e.,  for  any  x  €  fZ3, 
there  exists  c  €  C1  such  that  when  t  — *•  oo, 

i>{t-x)  =  c - - —  +  O(^).  (9) 

The  following  theorem  is  well-known,  and  is  a  direct  consequence  of  the  Gauss  theorem. 

Theorem  2.1 

Suppose  that  D  C  D  are  two  balls  in  R3,  and  that  D  is  bounded  by  a  sphere  S.  Suppose 
further  that  /  :  R3  \  D  — *  C  is  a  radiation  field  satisfying  the  equation  (6))  in  R3  \  D  and  the 
radiation  condition  (9))  at  oo.  Then  there  exist  two  analytical  functions  <7,77 :  S  — *■  C,  such 
that  for  all  x  €  R3  \  D , 

/(*)  =  jg  o’OO  *  fs(x)ds>  (10) 

and 

/(*)  =  fs  Vis)  •  fs,N{a)ix)ds,  (11) 

where  N (s)  denotes  the  exterior  normal  to  S  at  the  point  s. 


2.3.  Spherical  Bessel  and  Hankel  functions.  In  agreement  with  standard  practice, 
we  will  denote  by  jm  the  spherical  Bessel  function  of  the  first  kind  of  order  m,  and  by  hm,  the 
spherical  Hankel  function  of  order  m.  As  is  well  known  (see,  for  example  [18]),  jm  are  analytic 
on  the  whole  complex  plane  for  all  values  of  m,  while  hm  have  a  branch  cut  along  the  negative 
real  axis,  and  become  infinite  at  the  origin.  The  asymptotic  behaviour  of  the  functions  jm ,  hm 
for  large  m  is  given  by  the  formulae 


and 


/  v  en+h  .  zn+l 

lim  hm(z)  ■  ■=  — - —  =  1 

in— 00  ’  y/2  .  (2n  4- 1)" 
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(see  [1],  9.3.1,  9.3.2,  9.1.3).  For  large  z  and  fixed  m,  the  asymptotic  behavior  of  hm(z) 

is  given  by  the  formulae 


TTllC  T 

=  -§)-o(  e-w~). 


(14) 


L  /  \  i-(z---Z)  Im(ZK 

*  •  hm(z)  —  e  '  2  a>  =  0(- j^|2— ) 

when  z  — *■  oo,  as  long  as  Im(z)  >  0  (see  [1],  9.2.5,  9.2.7). 

All  spherical  Bessel  functions  are  ‘elementary  functions’.  In  particular, 


(15) 


ho(z)  =  - 


t  •  er 


1  _  *«*(*)  «<»(*) 

-  — - — 

*iW  =  + 

^jo(«)  =  -iiW, 

j 

M*)  =  -M*)* 


(16) 


dz 

The  following  theorem  is  known  as  the  Addition  Theorem  for  spherical  Bessel  functions,  and 
is  one  of  principal  analytical  tools  of  this  paper.  It  can  be  found,  for  example,  in  [1]. 


Theorem  2.2 

Suppose  that  r,  p,  8,  A  are  arbitrary  complex  numbers,  ij  =  ir  — 8,  and  that  R  €  C  is  defined 
by  the  formula 

R  =  (r2  +  p2  —  2  •  r  •  p  •  cos(0))»  =  (r2  +  p2  +  2  •  r  •  p  •  cos(r)))*  (17) 


(see  Figure  1).  Then 

£ (2»  +  1)  ‘  i«(A  •  r)  •  ;„(A  •  p)  •  Pn(cos(0))  = 

n=K) 

]T(-l)n  •  (2n  +  1)  •  j„(A  •  r)-jn( A  -p)  •  Pn(cos(»/))  = 


(18) 


n=0 
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(19) 


and 

E(2 n  +  1)  •  el*"'-'  •  jn(p)  •  P„(cos(0))  = 

n=0 

E(2n  +  1) .  in  .  Up)  .  Pn(cos(0))  = 

n=0 

If,  in  addition,  |  r  •  e±*'5  |<|  p  |,  then 


E(2n  +  1)  •  jn(A  •  r)  •  hn( A  •  p)  •  Pn(cos(0))  = 


n=0 

oo 


E(-l)n  •  (2n  +  1)  •  U A  •  r)  •  hn(A  •  p)  •  P„(co5(^))  = 


n=0 


ho(X  •  P)  = 


— *  •  e 


,i\R 


A  R 


(20) 


2.4.  Integrals  of  spherical  harmonics. 

A  function  u  :  S2  — *■  C  is  referred  to  as  a  spherical  harmonic  of  degree  n  if  the  function 
/  :  R?  — »  C  defined  by  the  formula 

f(x,  y,  z)  =  w(0,  <f> )  •  pm  (21) 

satisfies  the  Laplace  equation  in  R3  (see,  for  example,  [10]). 

Remark  2.1 

As  is  well-known  (see,  for  example,  [11]),  for  any  integer  n  >  0,  there  exist  exactly  2n  +  1 
linearly  independent  spherical  harmonics  of  order  n,  and  a  standard  representation  of  a  spherical 
harmonic  of  order  n  is  by  an  expression 

W(M)=  E  li-P?(cos(e))-eim  +,  (22) 

m=r— n 

with  P™  the  associated  Legendre  function  of  degree  n  and  order  m  (see,  for  example,  [1]),  and 
7 i,  *  =  0,  ±1,  ±2,  •  •  • ,  ±n  a  finite  sequence  of  complex  numbers.  However,  in  this  paper  we  will 
not  be  using  the  representation  (22),  remembering  only  that  the  spherical  harmonics  of  order 
n  constitute  a  complex  linear  space  of  dimension  2n  +  1. 
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We  will  need  the  following  three  well-known  lemmas  involving  the  integration  of  spherical 
harmonics  over  the  surface  of  the  sphere.  Lemmas  2.1,  2.2  below  can  be  found,  for  example,  in 
[10]),  [1],  respectively.  Lemma  2.3  is  a  simple  consequence  of  Lemma  2.2,  and  can  be  found  in 
[11- 

Lemma  2.1. 

For  any  spherical  harmonic  Y  of  degree  n  >  0, 

Y(s) '  Pn{cos(0{s)))ds  =  •Y(sn)  (23) 

Lemma  2.2. 

For  any  n  >  0  and  z  €  C, 

jn(z)  =  ei  *-",(®(4))  •  Pn(cos{0{s)))ds.  (24) 

The  following  theorem  is  a  simple  consequence  of  the  proceeding  two  lemmas,  Theorem  2.2, 
and  the  formulae  (12),  (13). 

Theorem  2.3. 

Suppose  that  p  6  R,  n  is  a  natural  number,  and  k  €  C  is  such  that  Im(k)  >  0.  Suppose 
further  that  u  €  R3  is  such  that  ||u||  <  p,  and  that  the  functions  Tn  :  S2  — ►  C,  Fn  :  R?  — ►  C 
are  defined  by  the  formulae 


n 


Tn(s)  =  Tn(0,  </>)  =  ]T  im  •  (2m  +  1)  •  hm(k  ■  p)  •  Pm{cos{0)), 

(25) 

m=0 

F„(tt)  =  - —  f  Tn(s)  ■  ei  k<E^ds. 

4  •  X  Js* 

(26) 

Then 

Fn(u)  =  M*  •  (p2  +  IMI2  +  2  •  p  •  INI  •  cos(i7))2 ), 

(27) 

with  Tf  the  angle  between  the  vector  u  and  the  z  axis.  Furthermore,  for  large  n, 

Proof. 


Combining  (25)  with  (26),  we  have 

P„(«)  =  •  f;  im  •  (2m  +  1)  •  hm(k  •  p)  .  f  Tn(s) .  Pm(cos(*))  •  eik<E^ds,  (29) 

4't  m=o  752 

and,  substituting  (19)  into  (29),  obtain 

Fn(u)  =  -i-.f^r*.(2m+l).M*-p)-  (30) 

4-7r  m=0 

[  Pm(cos(0(s)))  •  f)(2/+l).e1:§i.i1(Ar.||tt||)-PKc(P(ti),s)ys 

^S2  7^ 


1=0 


For  any  m  ^ 


f  Pm(cos(0(s)))  •  P,(c(P(u),s))ds  =  0,  (31) 

./S2 

since  in  this  case  Pm(cos(0(.s))),  Pj(c(P(tt),  s))  are  two  spherical  harmonics  of  different  degrees, 
and,  therefore,  orthogonal  on  S2.  Combining  (30)  with  (31),  we  have 

Pn(tt)  =  T^--i;(2m+l)2.Mfc-p)-im(*-|NI)  (32) 

4  ■  *  m=0 

•  /  J>m(cos(0(s)))  •  Pm(c(P(u),  s))ds. 

./S2 

Due  to  Lemma  2.1, 

jk  Pm(cos(<?(s)))  •  Pm(c(P(u),  s))ds  =  •  Pm(cos(0(P(u))),  (33) 

and  (32)  assumes  the  form 


F„(«)  =  2  (2m  +  !)  *  M*  •  *>)  *  im(*  •  IMI)  *  Pm(cos(0(P(u))).  (34) 

m=0 

Now  (27)  follows  from  the  combination  of  (34)  and  Theorem  2.2,  and  (28)  follows  from  the 
combination  of  (27)  with  (12),  (13). 


2.5.  Partial  wave  expansions  of  radiation  fields.  Suppose  that  the  function  4>  : 
R 3  — +  Cl  satisfies  the  Helmholtz  equation  (6)  outside  am  open  ball  D  of  radius  R  with  the 
center  at  the  point  Xo  €  i£3,  and  also  satisfies  the  radiation  condition  (9)  at  oo.  Then  there 
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exists  a  unique  sequence  of  spherical  harmonics  a  =  {am},  m  =  0, 1,2,  •  •  •,  such  that  for  any 
x  £  R3  \  D, 

<KX )  =  E  Q”*(5)  '  hm(kp),  (35) 

m—0 

with  ( p ,  s)  the  spherical  coordinates  of  the  vector  x—xq,  and  for  each  m  £  [0,  oo),  am  a  spherical 
harmonic  of  degree  m. 

If  a  function  ip  satisfies  the  equation  (6)  inside  D,  then  there  exists  a  unique  sequence 
of  spherical  harmonics  (3  =  {/?m}, m  =  0, 1,2, •••,  such  that  for  each  m,  /3m  is  a  harmonic  of 
degree  to,  and  for  any  x  £  D, 

V>(*)  =  53  A»(*)  3m(kp).  (36) 

m= 0 

A  derivation  of  the  formulae  (35),  (36)  can  be  found,  for  example,  in  [11],  and  we  will  refer 
to  functions  satisfying  the  Helmholtz  equation  as  radiation  fields,  to  expansions  of  forms  (35), 
(36)  as  h-expansions  and  j-expansions  respectively,  and  to  the  point  xo  as  the  center  of  the 
expansions  (35)  (36). 

The  following  lemma  is  a  direct  consequence  of  the  formulae  (12),  (13).  It  establishes  the 
convergence  rates  of  the  expansions  (35),  (36). 

Lemma  2.4 

If  D\  C  D  is  a  ball  of  radius  R\  <  R  with  the  center  at  xo  then  there  exists  c  >  0  such  that 
for  any  x  £  D\  and  N  >  |&|  •  Ri, 

W»)  -  E  <  (^)N-  (37) 

m=0  " 

If  D2  D  D  is  a  ball  of  radius  R2  >  R  with  the  center  at  xo  then  there  exists  c  >  0  such  that 
for  any  x  £  R?  \Dj  and  N  >  |fc|  •  R, 

\ip(x)  -  £;  (3m(0,  <p)hm(kp)\  <  (38) 

m= 0  * 

Remark  2.2 
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In  numerical  calculations,  expansions  (35),  (36)  are  truncated  after  a  finite  number  of 
terms,  and  the  resulting  expressions  are  viewed  as  approximations  to  the  fields  <f>,  rfr.  If  we 
want  to  approximate  4>  by  an  expansion  of  the  form  (37)  with  an  accuracy  e  then  according 
to  the  above  lemma,  we  have  to  choose 


,,,,  D  -  ln(e)  +  ln(c) . 

JV  > 

Since  logarithm  is  a  very  slowly  growing  function,  for  medium  and  large  scale  problems, 


(39) 


max(fi]  •  |fc|. 


—  ln(e)  +  ln(c) 
ln(J2)  —  ln(jRi ) 


(40) 


i.e.  the  number  of  terms  in  the  approximation  is  almost  independent  of  c,  and  must  be  roughly 
equal  to  |fc|  R\.  A  similar  calculation  shows  that  for  medium  to  large  scale  problems,  the 
expansion  (36)  can  be  truncated  after  approximately  N  >  \k\R  terms. 


2.6.  Numerical  integration  on  52. 


In  this  subsection  we  formulate  two  lemmas  describing  the  optimum  quadrature  formulas 
for  two  situations:  smooth  functions  on  a  circle,  and  smooth  functions  on  an  interval.  Then, 
we  use  these  lemmas  to  construct  a  high-order  quadrature  formula  on  S 2  (Theorem  2.4  below). 
Both  Lemmas  2.5,  2.6  are  well-known,  and  can  be  found,  for  example,  in  [15]. 

Lemma  2.5 


For  any  integer  m,  n  such  that  n  >  2|m|,  the  n-point  trapezoidal  quadrature  rule  on  the 
interval  [0, 2ir]  integrates  the  function  exmx  exactly. 

Lemma  2.6 


For  any  natural  n,  there  exist  a  unique  pair  of  finite  sequences  {x"),  {«>"},  *  =  1, 2,  •  •  •,  n, 
such  that  for  any  integer  k  6  [1, 2n  —  1], 

E»”(x?)‘=  /'«**•  (41) 

«=1  */-1 

Furthermore,  x"  €  [-1, 1]  and  w ?  e  [0, 1]  for  all  i  =  1, 2,  •  •  -n. 
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The  points  x"  and  the  coefficients  tu"  axe  known  as  the  nodes  and  coefficients  of  the  n-point 
Gaussian  quadrature  rule,  which  is  the  unique  n-point  quadrature  rule  that  integrates  exactly 
all  polynomials  of  order  up  to  2n  —  1. 

For  a  natural  n,  we  will  define  a  finite  sequence  “y^n  by  the  formula 

=  1).  (42) 

74 

and  the  finite  sequence  #i,02»  •  *  hy  the  formula 

dj  =  arccos(Wj).  (43) 

Now,  we  define  a  discretization  of  Dn  C  S2  as  a  collection  of  n2  points  s”fc  defined  by  the 
formula 
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3.1.  Sequences  of  spherical  harmonics  and  functions  on  52.  We  will  denote  by  Y 
the  set  of  all  sequences  a  =  {am},  m  =  0, 1,2,  •  •  •,  such  that  for  each  m,  am  is  a  spherical 
harmonic  of  degree  m.  We  will  define  a  norm  on  Y  via  the  formula 

ii«ii  =  («) 

m=0 

denote  by  X  the  subspace  of  Y  consisting  of  such  sequences  a  that  ||a||  <  oo>  and  observe  that 
the  norm  (48)  converts  X  into  a  Hilbert  space.  For  a  real  number  r  >  0,  we  will  denote  by  Xr 
the  subspace  of  X  consisting  of  all  sequences  a  =  {am},  m  =  0, 1, 2,  •  •  • ,  such  that 

Orn 

IK.II  •  (— r  •  -M  <  c  (49) 

er 

for  all  m  >  r.  We  will  denote  by  Yr  the  subspace  of  Y  consisting  of  all  complex  sequences 
/?  =  {/?m},  n*  =  0, 1, 2,  •  •  •,  such  that  for  some  c  >  0, 

IKUH^r-v/S^  (50) 

for  all  m  >  r.  It  is  easy  to  see  that  Xr  C  Yr,  and  that  the  condition  (49)  is  a  very  restrictive 
one,  since  in  order  to  satisfy  it,  the  elements  of  the  sequence  {am}  must  decay  roughly  as 
(r/2)m/m!,  while  the  condition  (49)  is  a  very  mild  one  -  it  prohibits  the  elements  of  {/3m}  from 
growing  faster  than  approximately  (2/r)m  •  m!.  By  applying  formulae  (9.3.1),  (9.3.2)  from  [1], 
it  is  easy  to  show  that  in  (35)  (36),  a  6  Yjfc|R  and  /?  €  Y\k\  r-  Conversely,  for  any  sequence 
a  €  Yw  r,  the  expansion  (35)  converges  inside  D,  and  for  any  /?  €  Vjfc|  R,  the  expansion 
(36),  converges  outside  D.  For  a  natural  n,  we  will  denote  by  Tn  a  linear  mapping  Yr  — *  Yr 
converting  a  sequence  a  =  {am},m  =  0, 1,2, •••  into  a  sequence  a  =  {dm},m  =  0,1,2,*”, 
defined  by  the  formulae 

«m  =  am  for  |m|  <  n 

am  =  0  for  |m|  >  n  + 1.  (51) 

Clearly,  Tn(Yr)  C  Xr,  and  for  obvious  reasons,  we  will  refer  to  T„  as  truncation. 

We  will  define  the  mapping  F  :  X  — *•  L2(S 2)  by  the  formula 

F(a)(s)  =  f)  am(s)  *  (52) 

m= 0 
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(53) 


with  a  =  {ao,ai, and  the  mapping  F_  :  X  — ►  £2(S2)  by  the  formula 
F-M 

m=0 

It  is  easy  to  see  that  the  mappings  F,  F_  axe  unitary  in  the  norm  on  X  defined  by  (48),  since 
the  expansion  into  spherical  harmonics  is  a  unitary  transformation,  and  any  two  spherical 
harmonics  of  different  degrees  are  orthogonal  to  each  other  (  see,  for  example,  [10])). 

The  following  obvious  lemma  can  be  found,  for  example,  in  [6].  It  connects  the  speed  of 
convergence  of  an  expansion  of  the  form  (52)  with  the  analyticity  of  its  sum. 

Lemma  3.1 

Suppose  that  a  G  XT  with  some  (arbitrarily  large)  r.  Then  F(a)  :  S2  — ►  C  is  an  analytic 
function  on  S 2. 

While  the  definitions  (52),  (53)  might  seem  arbitrary,  they  are  motivated  by  the  following 
two  lemmas,  which  are  a  direct  consequence  of  the  formulae  (14),  (15). 

Lemma  3.2 

If  <f> :  f£3  \  D  — »  C1  is  defined  by  (35),  then 

lim  <f>(x o  +  t  •  E(s))  •  t  •  e~t  k  t  =  F(a)(s).  (54) 

<—►00 


Lemma  3.3 

Suppose  that  :  D  — ►  Cl  is  defined  by  (36),  and  that,  in  addition,  (3  €  Xa  for  some 
(arbitrarily  large)  a.  Then 

lim  t  •  (tf(*o  +  t  •  E(s))  -  (F(j3)(s)  ■  ei  kt)  +  F_(/3)(s)  •  e~i  k  t).  (55) 

<—►00 


Remark  3.1 

The  above  two  lemmas  can  be  viewed  as  describing  the  far-field  behavior  of  the  potentials 
<f> ,  in  terms  of  the  mappings  F(a),  F((3)  :  52  -*  C,  and  we  will  refer  to  F(a),  F((3)  as  far-field 
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representations  of  <j>,  ij),  respectively.  Alternatively,  we  will  be  calling  F(o),  F(0)  far-field  forms 
of  expansions  (35),  (36). 


For  a  point  so  €  52,  a  natural  n,  and  a  complex  z,  we  will  define  the  function  A^n  :  S2  — *•  C1 
by  the  formula 

A^n(s)  =  5^  *m  •  (2m  +  1)  •  Pm(c(so,  s)  •  jm(z).  (56) 

m=0 

It  immediately  follows  from  (19)  that 

E  im  *  (2m  +  1)  .  Pj(c(s0,  s))  ■  jm(z )  =  (57) 

m=0 

and  for  »  =  oo,  (56)  assumes  the  form 

A^°°(s)  =  j  (58) 

For  a  point  so  €  52,  a  natural  n,  and  a  complex  z ,  we  will  define  the  function  :  S2  — *•  C1 

by  the  formula 

A‘«”(5)  =  E  *w  ’  (2m  +  !)  ‘  ^m(c(so,s))  •  hm(z),  (59) 

m=0 


and  observe  that  no  analogue  of  the  formula  (58)  is  possible  in  this  case  (at  least  in  the  proper 
sense),  since  the  series 

E  *'m  •  (2m  +  X)  •  PM* 0,  s))  •  hm(z)  (60) 

m=0 

does  not  converge. 

Finally,  we  will  define  mappings  A^n,  Af£n  :  L2(S2)  — ►  L2(S2)  via  the  formulae 

as"(/)« = ■*:;»  •/(>).  (si) 

"£"(/)«  =  rfTM/M  (62) 

respectively,  with  /  €  L2(S2). 

3.3.  Definition  of  translation  operators.  For  the  remainder  of  this  section,  D\,  D?, 
Dz  will  denote  three  balls  in  R3  such  that  D2  C  D\  and  D\  O  D3  =  0  (see  Figure  2).  The 
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centers  and  radii  of  these  disks  will  be  denoted  by  ei,  C2>  C3  and  R\,  R2,  R3  respectively.  We 
will  denote  the  spherical  coordinates  of  the  vector  C2  —  cj  by  (pu,  ^12),  the  spherical  coordinates 
of  the  vector  Ci  —  C2  by  (p2i»S2i)»  the  spherical  coordinates  of  the  vector  C3  —  Cj  by  (P13,  £13), 
and  the  spherical  coordinates  of  the  vector  cj  —  C3  by  (/»3i,S3i).  For  a  point  x  G  f£3,  we  will 
denote  by  (pi,$i),  (p2,$2),  (P3, S3)  its  spherical  coordinates  with  respect  to  the  centers  c\,  C2, 
C3  respectively. 

Suppose  now  that  xj) :  R3  -+  C1  is  a  radiation  field  analytical  in  R?  \  D 2  and  satisfying  the 
radiation  condition  (9)  at  00.  Suppose  further  that  xj}  is  represented  by  an  expansion 

^(x)  =  5Z  2)  •  hm(kp2)  (63) 

m= 0 

valid  in  i?3  \  Z?2>  and  by  an  expansion 

V’(z)  =  S  &n(5l)  •  hm(lepi)  (64) 

m=0 

valid  in  R3\D\.  It  is  easy  to  see  that  (3  G  A|jt|.R,  depends  linearly  on  (3  G  and  we  will 

denote  by  UC2,c\  the  operator  such  that 

P  =  Ue 2<ei(P).  (65) 

Suppose  that  <f> :  R3  — ►  C1  is  a  radiation  field  analytical  in  D\  and  represented  by  an  expansion 

<K*)=jt°m(s1).jm(kp1)  (66) 

m=0 

valid  in  D\,  and  by  an  expansion 

<KX )  =  &xn(s2)  • jm(kp2 ),  (67) 

m=0 

valid  in  Z>2-  Again,  it  is  easy  to  see  that  a  G  Fj*|.R2  depends  linearly  on  a  G  Yjfc|.R,,  and  we 
will  denote  by  Vc\fi2  the  operator  Fjk|.Rj  —*  yjfc|.R2  such  that 

6  =  Vcl,c2(a).  (68) 

For  any  r  >  0,  we  will  denote  by  Kl.c2  the  restriction  of  VcitC2  on  the  subspace  Xr  of  , 
so  that 

=  (VcU1)IXr.  (69) 
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Finally,  suppose  that  ip  :  R3  \  D\  — ►  C1  is  a  radiation  field  analytical  outside  the  ball  D\  and 
satisfying  the  radiation  condition  (9)  and  that  it  is  represented  in  R 3  \  D\  by  the  expansion 
(63).  Then  inside  the  ball  £>3,  the  function  ip  can  be  represented  in  the  form 

^(*)  =  E  7m(«8)  •  Jm(*As)  (70) 

m= 0 

with  7  €  yjjc|./j3  a  linear  function  of  a  €  and  we  define  the  operator  We xtC2  ■  — *• 

yjfc|.R3  via  the  formula 

7  =  W*,c3(a).  (71) 

3.4.  Diagonal  Forms  of  Translation  Operators.  This  subsection  describes  the  diagonal 
forms  of  the  translation  operators  U,  V,  W  for  the  Helmholtz  equation.  These  diagonal  forms 
are  provided  by  the  Theorems  3.1  -  3.3  below,  and  are  the  principal  purpose  of  this  paper. 

Theorem  3.1. 

If  the  operator  UC2,ci  '•  -X)*) .r3  — ►  is  defined  by  the  formula  (65),  then 

UC2,cl  =  F~l  o  A*'?12’00  o  F.  (72) 


Proof. 

We  will  prove  (72)  by  showing  that 

FoHd|Cl=Aj;fwoF.  (73) 

Suppose  that  s  €  S2,  and  (3  6  Combining  (61)  with  (58),  we  have 

A^*’00(s)  =  F(p(s)  •  ei  k  R»-e^-^s\  (74) 

On  the  other  hand,  due  to  Lemma  3.1, 

F0)(s)  =  lim  +  t  •  £(«))  •  *  •  e~i  k  t 

t — ►OO 

=  lim  <P(C2  +  (C!  -  c2)  +  t  ■  E(s))  ■  t .  e~i  k  t  (75) 

t—+ao 

=  fc  + ]j«+7^i  •  u« + <  •  »n)  •  <  • 
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(76) 


with  u  =  ci  —  C2,  and  v  =  E(s).  Denoting  ||u  +  t  •  r||  by  r,  we  obtain  after  simple  analysis 
t  =  r  -  («,w)  +  C>(-), 


and  (75)  assumes  the  form 


F(0)(s)  =  lim#C2  +  -  —  ~  — — — ^(t))  w  .r) 
(  iK  }  ™ =°n  +  Hu +  (r  _(«,«)  + 0(1)).  ®||  j 

.  (r  -  («,  v)  +  0(1))  •  e-»-fc-(T-(u,v)+0(i )) } 


(77) 


which,  due  to  the  radiation  condition  (9)  can  be  reduced  to 

r-oo  ||u  + (r- (u,v))-v|| 

•  r  •  ei‘k'(u'v) .  e~i  k'T. 


(78) 


Finally,  combining  Lemma  3.1  with  (78)  yelds 

F0)(s)  =  Km  ^c2  +  |j^.r).r.e-<^.e*fe(M 

=  lim  #c2  +  E(s)  ■  t)  ■  t  •  e— *  T  •  e*^*1  (79) 

T  — ►OO 

=  F(/3)(s)  •  e*'fc‘(c,-C2,£'(*)). 

Now,  the  conclusion  of  the  theorem  follows  from  the  combination  of  (79)  with  (74). 

The  above  theorem  provides  the  diagonal  form  of  the  operator  Uc^,Cl  shifting  the  origin 
of  an  h  -  expansion.  Theorem  3.2  below  is  the  analogue  of  of  Theorem  3.1  for  the  case  of 
j-expansions.  Since  the  proofs  of  the  two  theorems  are  virtually  adentical,  we  omit  the  proof 
of  the  following  theorem. 

Theorem  3.2. 

If  the  operator  VcitC2  :  -*•  Fj*|..Rj  is  defined  by  the  formula  (68),  then  for  any  r  >  0, 

Vci,c2  =  F-1oAk'f»’°°oF.  (80) 


The  following  two  lemmas  are  an  immediate  consequence  of  Theorems  3.1,  3.2.  Lemma  3.4 
provides  the  far-field  representations  of  the  potentials  (7),  (8)  of  charge  and  dipole.  Given  a 
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far-field  representation  of  a  potential  of  the  form  (36),  Lemma  3.5  provides  an  expression  of 
its  value  (and  the  value  of  its  gradient)  at  any  point  within  the  region  of  its  validity. 

Lemma  3.4 


Suppose  that  in  (35), 


4>  =  4>kXi, 

with  x\  an  arbitrary  point  in  R3.  Then  for  any  s  £  S2, 
F(ox>) = 

If  h  €  .ft3  is  such  that  ||/t||  =  1,  and 
<t>  =  4>xuhi 


then  for  any  $  6  52, 

F(a)(s) 


aH|xo-x1||,°0(5)  . . .  (so-gi.ffi) 

Pixo-xt)  v  )  11*0  -*l|| 

eiMxo-*i ,£(.))  .i.k.  c(p(X()  _  Xl)>s). 


(81) 


(82) 


(83) 


(84) 


Lemma  3.5 


Suppose  that  the  potential  ip  is  defined  by  (36),  and  that  /?  €  Xr,  with  some  0  <  r  <  oo. 
Then  for  any  x  e  D  and  h  6  R3  such  that  ||h||  =  1, 

iP(x)  =  (  F(/3)(s)  ■  A pfxoZx)'00^  =  [  F(P)(*)  ■  eik<*°-*’EWds.  (85) 

and 


d_ 

dt 


V>(*  +  t  •  h)\t=0  =  f  F(0)(s)  •  Ap{l*“_*)l|,0O(s)  -i  k 

J  s* 


(*o  -  x ,  £(s)) 

11*0  -  *11 


ds 


=  /  F(/3)(s)  •  ei  k  (*o-x£(*)) .  i .  Jt .  c(P(x0  -  *),  s)ds. 
Js* 


(86) 


While  the  preceeding  two  theorems  are  fairly  obvious,  and  appear  to  be  known  (though 
in  a  somewhat  different  form)  among  certain  groups  of  physicists,  the  following  theorem  is 
considerably  more  technical,  and  appears  to  be  quite  new. 
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Theorem  3.3. 


Suppose  that  the  operator  Wci,c3  •  T|fc|.rt3  is  defined  by  the  formula  (71).  Suppose 

further  that  rp  :  R3  \  D\  — *  C  is  a  radiation  field  represented  by  the  expansion  (64)  outside 
Di,  and  by  the  expansion  (70)  inside  D3.  For  any  n  >  1,  we  will  denote  by  ipn  the  radiation 
field  D3  — ►  C  defined  by  the  formula 

V>n(*)  =  £  7m(*»)  •  jm(kf>3),  (87) 

m= 0 

with  7n  =  {7",  7J ,  73 ,  •  •  •}  defined  by  the  formula 

7n  =  F~l  o  Jf  *; f3-n  o  F(a).  (88) 

Then  for  any  x  €  D3, 

lim  rpn(x)  =  ip(x).  (89) 

n— *00 

Furthermore, 

max  |  V>„(x)  -  #(x)  |=  0((^^)"  -  Ml).  (90) 

Proof. 

Due  to  Theorem  2.1,  it  is  sufficient  to  prove  (90)  in  the  special  case  of 

=  fx0,  (91) 

with  Xq  an  arbitrary  point  in  D\.  Combining  (91)  with  (82),  we  have 

F((3)(s)  =  €ik<Xo~Ci  'EW  (92) 

for  all  s  €  S2,  and  combining  (92)  with  (62),  (59),  obtain 

M«?'n  °  rm*)  =  (93) 

eik(x0-cuE(a)) .  £  im  ■  (2m  +  1)  •  Pm(c(sl3,s ))  •  hm(k  ■  pl3).  (94) 

771=0 
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Now,  combining  (93)  with  (85),  (87),  we  have 


V’n(x)  = 


ei-k-(xo-ci,E(s))  m  ei  k  (x-c3,E(s))  m 


£  ,m  '  (2m  +  !)  •  pm(c(si3,  s))  •  hm(k  ■  P13), 

m=0 


(95) 


and  (90)  follows  from  the  combination  of  (95)  and  Theorem  2.3. 

3.5.  Numerical  evaluation  of  translation  operators.  For  the  rest  of  this  paper, 
we  will  view  the  asymptotic  representations  F( or),  /(/?)  defined  by  (52)  (as  opposed  to  the 
expansions  of  the  forms  (35),  (36))  as  our  principal  tool  for  representing  radiation  fields. 
Lemma  3.4  permits  one  to  calculate  asymptotic  representations  of  fields  of  distributions  of 
charges  and  dipoles  without  evaluating  the  coefficients  of  their  h-expansions,  and  Lemma  3.5 
provides  a  tool  for  calculating  the  fields  and  derivatives  of  the  fields  with  given  asymptotic 
representations  without  having  to  evaluate  the  coefficients  of  j -expansions  of  these  fields. 

For  a  radiation  field  if) :  R3  —+  C1  analytical  outside  D\  given  by  the  expansion  (64),  and 
an  integer  n  >  2,  we  will  denote  by  F^Ci  the  function  F0)  tabulated  at  the  n2  nodes  s"fc 
defined  by  the  formulae  (42)  -  (44),  so  that 


,k)=Fmsi„). 


(96) 


Similarly,  for  a  radiation  field  <f>  analytical  inside  D\  defined  by  (66)  and  possessing  an  asymp¬ 
totic  representation  F(a),  we  will  denote  by  G1^  Cj  the  table  of  n2  complex  numbers  defined  by 
the  formula 


Gn^Cl(j,k)  =  F(a)(slk).  (97) 

and  view  ,  G^  CJ  as  finite-dimensional  projections  of  the  asymptotic  representations  of  the 
radiation  fields  Vs  4>- 

Given  a  function  G  :  Dn  — ►  C  (see  (44),  we  will  consider  a  radiation  field  G  :  R3  — »  C1 
defined  by  the  formula 

G(*)=  E  *)■£-*•<*•-■*<•«)>.  (98) 

j,k=l 
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Clearly,  (98)  is  a  quadrature  formula  approximating  the  integral  (85)  and  we  will  look  upon 
(98)  as  an  approximation  to  the  field  <f>.  Differentiating  (98)  with  respect  to  x,  we  obtain  the 
formula 

jtG(i  +  <*)(,_*  = 

<  •  *  •  IM  •  X>J*  •  GJ,tl  0,  *)  •  c(P(x„  -*),.)•  (99) 

j,k=l 

with  any  h  €  R3-  Finally,  we  will  define  mappings  S™*  :  Cnxn  — ►  CnXT*  by  the 

formulae 


pmn  _  a  fc-pi2,m 

xC2Ci  i  C2,Ci 


(100) 


£)mn  _  a  fc*pi2,m 
^CiC2  ~~  iXC\,C2 


(101) 


^cic3  iWci,c3  |j 


(102) 


Observing  that  ,  S™”  are  restrictions  to  the  nodes  of  the  discretization  (45)  of 

the  mappings  A*'^’”*,  Af*'^3,m  respectively,  we  will  look  upon  the  operators  , 

S™”3  as  approximate  discretizations  of  diagonal  forms  of  the  operators  £/££ ,  Fc™",  W7™" . 

Remark  3.5 

By  combining  the  above  lemma  with  Remark  2.2,  it  is  easy  see  that  the  number  n  of  nodes 
in  the  discretization  GJ  of  the  function  G^,Cl  :  S2  — ►  C1  has  to  be  approximately  equal 
to  (2  •  |fc|  •  f?j)2,  and  is  almost  independent  of  the  accuracy  e  with  which  the  field  4>  is  being 
calculated. 

Theorems  3.1  -  3.3  provide  a  tool  for  shifting  the  origins  of  asymptotic  expansions  of  radi¬ 
ation  fields,  and  for  converting  asymptotic  representations  of  the  form  (54)  into  asymptotic 
representations  of  the  form  (55)  for  a  cost  proportional  to  n,  where  n  is  the  number  of  nodes 
in  the  discretization  (45)  of  the  interval  S2.  In  the  following  two  sections,  this  apparatus  is 
used  to  construct  an  algorithm  for  rapid  evaluation  of  radiation  fields  of  charge  (and  dipole) 
distributions. 


21 


4.  Rapid  Evaluation  of  Radiation  Fields  of  Charge  Distributions 

In  this  section,  we  describe  an  algorithm  for  rapid  evaluation  of  the  field  and  the  normal 
derivative  of  the  field  created  on  a  surface  T  C  -R3  by  charge  and  dipole  distributions  on  that 
same  surface.  For  definitiveness,  we  will  be  discussing  the  evaluation  of  the  field  created  by 
a  charge  distribution.  The  algorithms  evaluating  the  normal  derivative  of  the  field  created  by 
a  charge  distribution,  and  the  field  and  the  normal  derivative  of  the  field  created  by  a  dipole 
distribution  are  quite  similar. 

4.1.  Notation.  We  will  consider  the  situation  depicted  in  Figure  3.  The  surface  r  C  R3 
is  discretized  into  n  nodes  xj,  X2,  •  •  •,  xn,  and  we  will  assume  that  these  nodes  are  distributed 
on  the  surface  in  a  roughly  uniform  manner.  Suppose  that  for  each  i  =  1,2,  ••*,«,  a  charge 
o,  of  strength  a,-  is  located  at  the  point  Xj.  In  this  section,  we  describe  an  algorithm  for  rapid 
calculation  of  approximations  <7,  =  1, 2,  •  •  ■ ,  n  to  the  sums 

G.(*  .')=  t,  "><(*■)  I103) 

i-i, 

for  i  =  1,2,  •••,».  Clearly,  this  is  an  order  n2  process  (evaluating  n  fields  at  n  points). 
However,  if  we  are  interested  in  evaluating  (103)  with  a  finite  accuracy  (which  is  always  the 
case  in  practical  calculations),  Theorem  3.3  and  Lemmas  3.4,  3.5  can  be  used  to  speed  up  the 
process. 

For  an  integer  m  >  4,  we  will  subdivide  the  surface  T  into  m  non-intersecting  patches  P,, 
each  patch  roughly  rectangular  in  shape,  and  containing  approximately  ^  of  the  nodes  x,. 

For  each  i  €  [l,m],  we  will  denote  by  r,  the  radius  of  the  smallest  circle  containing  P, ,  and 
define  r  6  R  by  the  formula 

t—  max  r,-.  (104) 

Remark  4.1 


Obviously,  for  sufficiently  large  m, 
L 
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(105) 


where  L  is  the  diameter  of  the  surface  T.  This  will  be  important  in  Subsection  4.3. 

For  each  t  =  1, 2,  •  •  • ,  m,  we  will  denote  by  Zj  the  center  of  the  i  —  th  patch,  by  Aj  the  set 
of  all  charges  aj  such  that  aj  E  A,-,  and  by  D:  the  disk  of  radius  r  with  the  center  at  Zj. 

We  will  denote  by  Wj  the  union  of  all  A,-  such  that  \\z:  —  z,  ||  >  3  •  r,  and  by  Wj  the  union  of 
all  Aj  such  that  || zj  —  z,-||  <  3  •  r.  Obviously,  Aj  C  Dj  for  any  j  =  1, 2,  •  •  • ,  m.  Also,  it  follows 
from  the  triangle  inequality  that 

min  ||x-y||>r  (106) 

j,  y€Aj 

for  any  i,j  such  that  Aj  C  Wj.  Finally,  we  will  denote  by  4>j  the  field  of  all  charges  Oj  such 
that  Xj  C  Aj  and  observe  that  if  xp  G  Aj  then 

Gcr(xp)=  Yi  <f>i(xP)+  Y  VitxiiZp)-  (107) 

AiCWj  xi&Wj 

4.2.  Detailed  description  of  an  order  jV3/2  algorithm.  In  this  subsection,  M,  N  will 
denote  "sufficiently  large”  integer  numbers.  The  actual  choice  of  the  numbers  M,  N  is  discussed 
in  the  following  subsection. 

We  will  evaluate  the  fields  (103)  in  five  steps. 

Step  1. 

Using  Lemma  3.4,  obtain  discretized  asymptotic  representations  of  the  fields  4>j  for 

all  j  =  1,2 
Step  2 

For  every  pair  of  natural  numbers  i,j  €  [1,  m]  such  that  Aj  C  Wj,  calculate  the  representa¬ 
tion 


of  the  field  V\j  =  6-j,  and  view  it  as  a  finite  -  dimensional  approximation  to  the  asymptotic 
representation  of  the  field  <f>i  on  Dj. 

Step  3 
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For  each  natural  j  €  [1,  m],  calculate  the  sum 


G 


N  _ 


A,CW} 


r>N 


(109) 


and  view  the  field  V’j  =  35  an  approximation  to  the  field  Y^AiCW,  Qn  and  G^>z.  as  a 

finite-dimensional  approximation  to  the  asymptotic  representation  of  ipj  on  Dj. 

Step  4 

For  each  natural  j  €  [1,  m],  evaluate 


(110) 


for  all  i  such  that  x,-  €  Pj  and  look  upon  (110)  as  an  approximation  to 
Step  5. 

For  each  j  =  1,2 ,  •  •  • ,  to,  evaluate  the  sum 


]C  °p^P(xi)  (HI) 

xpe  Wj 

for  all  i  such  that  x,-  £  Pj,  and  view  (111)  as  an  approximation  to  G<,(xi). 

4.3.  Choice  of  parameters  and  CPU  time  estimate.  In  the  estimates  below,  a,  b,  c,  d,  e 
are  coefficients  determined  by  the  computer  system,  language,  particular  implementation  of  the 
algorithm,  etc. 

Step  1 

Obviously,  this  step  will  require  order  n  •  Ar2  operations  (tabulating  F^.  ^  at  N2  nodes  on 
Dn  for  each  of  the  nodes  xi,X2, •••,!„  ).  Combining  Remarks  3.5,  4.1,  we  observe  that  N  ~ 
\k\-L/y/m ,  and  the  CPU  time  estimate  for  this  step  becomes  a-n-{\k\-L/m)2  =  a-n-{\k\-L)2 [m. 

Step  2 

For  each  of  the  pairs  i,j  such  that  Aj  C  Wi,  evaluating  (108)  will  require  order  N2 
operations  (see  (102)  ),  and  the  total  number  of  such  pairs  is  less  than  m2,  which  results  in 
the  CPU  time  estimate  of  b  •  m2  •  n  ~  b  •  m2  •  (|fc|  •  L/y/m)2  =  b  ■  m  •  (| Ar|  •  L)2  for  this  step. 

Step  3 

Obviously,  evaluating  the  sums  (109)  for  all  j  =  1, 2,  •  •  • ,  m  is  an  order  c  ■  m  ■  N2  ~ 
c  •  m  •  (|fc|  •  L/y/m)2  =  c  •  (|fc|  •  L)2  procedure. 
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Step  4 

Evaluating  (110)  for  each  *  =  1,2,  •  •  -  ,n  is  an  order  N2  procedure,  resulting  in  the  total 
CPU  time  estimate  for  this  step  of  d  •  n  •  N2  ~  d  •  n  •  (|*|  •  L)2/m. 

Step  5 

Evaluating  the  sum  (111)  for  each  t  =  1,2,  is  an  order  n/m  procedure,  with  the 
resulting  CPU  time  estimate  of  e  •  n2/m  for  this  step. 

Summing  up  the  time  estimates  for  the  steps  1-5,  we  obtain  the  following  time  estimate  for 
the  whole  process: 


T  =  — n'(|fc|'X)2  +  b  •  m  •  (| *|  •  Lf  +  c  •  (|Jb|  •  Lf  + 


€  •  71 


(112) 

m  m 

with  A  =  a+d ,  and  we  would  like  to  choose  m  in  such  a  manner  that  (112)  would  be  minimized. 
Differentiating  (112)  with  respect  to  m,  and  setting  the  resulting  derivative  to  zero,  we  obtain 


/A  •  n  •  (|A;|  •  X)2  +  e  •  n2 

mmin  =  V - Ml*i~:i)* - 


(113) 


and  the  corresponding  minimum  of  (112)  is  equal  to 


Tm m  =  2  y/A  •  n  .  (|fc|  L)2  +  e  •  n2  •  ^6  •  (|*|  •  L)2  +  c  •  (|*|  •  Lf.  (114) 


If  the  calculations  are  performed  with  a  ^xed  number  of  nodes  per  wavelength  (which  is  often 
a  reasonable  assumption),  n  is  proportional  to  (|*|Z)2,  and  (114)  assumes  the  form 


Tm ^ 


(1*1 


(115) 


or 

Tjpjn  ~  n*,  (116) 

which  for  large  n  is  considerably  smaller  than  n2. 

4.4.  Further  reduction  of  the  CPU  time  estimate  of  the  process.  The  approach 
of  the  above  subsection  can  be  used  recursively  by  subdividing  each  of  the  sets  A i  into  subsets 
=  1,2,  with  appropriately  chosen  m  and  representing  the  fields  fa  as  sums  fa  = 

5Z,  fa j  where  faj  is  the  field  created  by  all  charges  ap  such  that  ap  €  Bij.  A  calculation  similar 
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to  the  one  in  the  preceeding  section  shows  that  such  an  algorithm  will  have  an  asymptotic  CPU 
time  estimate  of  n4/3. 

By  continuing  this  process  recursively  until  only  a  finite  number  of  nodes  is  left  on  a  surface 
patch  on  the  finest  lavel,  one  can  obtain  an  order  nlog(n)  algorithm  for  evaluating  (103). 
However,  our  estimates  indicate  that  for  problems  of  practicable  size  (n  <  1000,000),  the 
improvement  in  actual  computation  times  obtained  by  replacing  an  order  n4/3  algorithm  with 
an  order  nlogn  algorithm  would  not  be  very  significant. 
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